Cache Oblivious Matrix Operations Using Peano Curves

Algorithms are called cache oblivious, if they are designed to benefit from any kind of cache hierarchy—regardless of its size or number of cache levels. In linear algebra computations, block recursive approaches are a common approach that, by construction, lead to inherently local data access pattern, and thus to an overall good cache performance[3]. In this article, we present block recursive approaches that use an element ordering based on a Peano space filling curve to store the matrix elements. We present algorithms for matrix multiplication and LU decomposition, which are able to minimize the number of cache misses on any cache level. 1 A block recursive scheme for matrix multiplication Consider the multiplication of two 3× 3-matrices, such as given in equation (1), where the indices of the matrix elements indicate the order in which the elements are stored in memory. a0 a5 a6 a1 a4 a7 a2 a3 a8  } {{ } =: A  b0 b5 b6 b1 b4 b7 b2 b3 b8  } {{ } =: B =  c0 c5 c6 c1 c4 c7 c2 c3 c8  . } {{ } =: C (1) The scheme is similar to a column-major ordering, however, the order of the even-numbered columns have been inverted, which leads to a meandering scheme, which is also equivalent to the basic pattern of a Peano space filling curve. Now, if we examine the operations to compute the elements cr of the result matrix, we note that the operations can be executed in a very convenient order – from each operation to the next, an element is either reused or one of its direct neighbours in memory is accessed: